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The material derivative concept of continuum mechanics and an adjoint 
variable method of design sensitivity analysis are used to relate 
variations in structural shape to measures of structural performance. A 
domain method of shape design sensitivity analysis is used to best utilize 
the basic character of the finite element method that gives accurate 
information not on the boundary but in the domain. Implementation of 
shape design sensitivity analysis using finite element computer codes is 
discussed. Recent numerical results are used to demonstrate accuracy that 
can be obtained using the method. Result of design sensitivity analysis 
is used to carry out design optimization of a built-up structure. 


If 2 1. INTRODUCTION 

M *-> I 

cn B- tvj 

cn m cn 

^ tj A substantial literature has been developed in the field of shape 

O 

design sensitivity analysis and optimization of structural components 
[1-3] over the past few years. Contributions to this field have been made 
using two fundamentally different approaches to structural modeling and 
analysis. The first approach uses a discretized structural model, based 
on finite element analysis, and proceeds to carry out shape design 
sensitivity analysis by controlling finite element node movement and 
differentiating the algebraic finite element equations [4-6]. The second 
approach to shape design sensitivity analysis uses an elasticity model of 




*Research supported by NASA-Langley Grant NAG-1-215. 


the structure and the material derivative method of continuum mechanics to 
account for changes in shape of the structure [7-13], Using this 
approach, expressions for design sensitivity in terms of domain shape 
change are derived in the continuous setting and evaluated using any 
available method of structural analysis; e.g., finite element analysis, 
boundary element analysis, photoelasticity, etc. 

Shape design sensitivity analysis for several structural components 
has been treated in Refs. 2, 9, and 10 where sensitivity information is 
explicitly expressed as integrals, using integration by parts and boundary 
and/or interface conditions to obtain identities for transformation of 
domain integrals to boundary integrals. Numerical calculation of design 
sensitivity information in terms of the resulting boundary integrals thus 
requires stresses, strains, and/or normal derivatives of state and adjoint 
variables on the boundary. However, when the finite element method is 
used for analysis of built-up structures, the accuracy of numerical 
results for state and adjoint variables on interface boundaries may not be 
good [14], 

To overcome this difficulty, a domain method of shape design 
sensitivity analysis is developed in Ref. 15, in which design sensitivity 
information is expressed as domain integrals, instead of boundary 
integrals (boundary method). The domain and the boundary methods are 
analytically equivalent. However, when one uses an approximate numerical 
method such as finite element analysis, the resulting design sensitivity 
approximations may give quite different numerical values. Moreover, the 
domain method offers a remarkable simplification in derivation of shape 
design sensitivity formulas for built-up structures since interface 
conditions are not required to obtain shape design sensitivity formulas. 

In the domain method, numerical evaluation of the sensitivity information 
is more complicated and inefficient than the result of the boundary 
method, since the domain method requires integration over the entire 
domain, whereas the boundary method requires integration over only the 
variable boundary. To alleviate this problem, a boundary layer of finite 
elements that vary during the perturbation of the shape of a structural 
component is introduced in Ref. 16. 

In shape design problems, nodal points of the finite element model 
move as shape changes. In Ref. 17, a method of automatic regridding to 
account for shape change has been developed using a velocity field in the 


domain that obeys the governing deformation equations of the elastic 
solid. 

Using the domain method of Ref. 15 and results of conventional 
(sizing) design sensitivity analysis theory of Ref. 2, a design component 
method is developed in Ref. 18 for unified and systematic organization of 
design sensitivity analysis of built-up structures, with both conventional 
and shape design variables. That is, conventional and shape design 
sensitivity formulas for each standard component type can be derived. The 
result is standard formulas that can be used for design sensitivity 
analysis of built-up structures, by simply adding contributions from each 
component. The method gives a systematic organization of computations for 
design sensitivity analysis that is similar to the way in which 
computations are organized within a finite element code. 

A numerical method has been developed in Ref. 19 to implement the 
results of the design component method, using the versatility and 
convenience of existing finite element codes. It is shown in Ref. 19 that 
calculations can be carried out outside existing finite element codes, 
using postprocessing data only. Thus, design sensitivity analysis 
software does not have to be imbedded in an existing finite element code. 

The purpose of this paper is to combine these developments to present 
a unified method of shape design sensitivity analysis and numerical 
implementation of the method with existing finite element codes. Even 
though only static response is considered here, the method is also 
applicable for eigenvalue design sensitivity analysis as shown in Ref. 2. 

The design sensitivity analysis method presented here supports 
optimality criteria method for structural optimization and serves as the 
foundation for iterative methods of structural optimization using 
nonlinear programming. At a more practical level, the design sensitivity 
analysis method can be used to develop an interactive computer-aided 
design system [2], A large scale built-up structure is optimized to 
demonstrate capability of the method. 


2. VARIATIONAL FORM OF GOVERNING EQUATIONS 

While a substantial library of structural components must be 
considered to implement the design sensitivity analysis for broad classes 
of applications, the component library to be considered in this paper is 


limited to truss, beam, plate, plane elastic solid, and three dimensional 
solid components. Even though this is a somewhat restricted class of 
components, it is general enough that significant applications can be made 
and practicality of the method can be demonstrated. 

In the actual formulation, the truss and beam components, including 
both bending and torsion of the beam, are incorporated into a single 
component. Similarly, plate and plane elastic solid components are 
combined as a single component. To be more specific, the following 
formation of beam/truss, three dimensional elastic solid, and plate/plane 
elastic solid components are employed: 

A. BEAM/TRUSS 

Consider the beam/truss component of Fig. 2.1. The energy bilinear 
form (internal virtual work) [2] of the component is 



i __ 

+ / hEv v dx [2.1] 

0 x * 
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Figure 2.1 Beam/Truss Component 

1 2 

where w , w , 0, and v are two orthogonal lateral displacements, angle of 

12 T 

twist, and axial displacement, respectively, and z = [w , w , 6, v] . 
Throughout this paper, an overbar; e.g., 7 , denotes a virtual displace- 
ment. Subscript x in Eq. 2.1 denotes derivative with respect to x. In 
Eq. 2.1, E, G, I*, I^, J, and h are Young's modulus, shear modulus, two 
moments of inertia, torsional moment of inertia, and cross-sectional area 
of the component, respectively. The conventional design variable is u = 
h(x) and the shape design variable is the length of the domain fl = [0,i], 
The load linear form (external virtual work) [2] of the component is 



[ 2 . 2 ] 


A* ■* 1 ** O O ** ** 

l (7) = / q w dx + / q w dx + / r0dx + / f7dx 

u ‘“ 0 0 0 0 

where q 1 , q 2 , r, and f are two orthogonal lateral loads, twisting moment, 
and axial load, respectively, as shown in Fig. 2.2 [2]. If there are 

point loads, a Dirac delta measure can be used for q*, q 2 , r, and f in 

Eq. 2.2 [2]. 


q 2 U) 



Figure 2.2 External Loads For Beam/Truss 

The variational equation of the beam/truss component is [2] 

a u>J} (z,7)= * u>n (z), for al 1 7 £ Z [2.3] 

where Z is the space of kinematically admissible displacement. That 
is, Zc [H 2 (0,a)] 2 x [H 1 (0,&)] 2 and elements of Z satisfy kinematic 
boundary conditions where H^O,*) is the Sobolev space of order i [2]. As 
possible boundary conditions, the beam/truss component can be simply 
supported, clamped, cantilevered, or clamped-simply supported. It is 
shown in Ref. 2 that the variational Eq. 2.3 is applicable for all 
boundary conditions mentioned. 


B. THREE DIMENSIONAL ELASTIC SOLID 

Consider the three dimensional elastic solid of Fig. 2.3. For plane 
elastic solid, results of the three dimensional elastic solid may be 
reduced. 

The strain tensor is defined as 


e iJ (z) -i (z] + z\), ij = 1,2,3, 


xg n 


[2.4] 


r 



Figure 2.3 Three Dimensional Elastic Solid 
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where z s [z , z , z ] is displacement field and subscript i, i = 1,2,3, 
denotes derivatives with respect to variable x^. The stress-strain 
relation (generalized Hooke's Law) is 


y lj (z) = l C ljk£ e kA (z), 1,j,k,i = 1,2,3, xe« 

k ,Jt=l 


C2.5] 


where C is the elastic modulus tensor, satisfying symmetry relations 
c ijk£ = c jik£ and c ijk£ = i ,j ,k , £ = 1,2,3. The energy bilinear 

form [2] of the three dimensional elastic solid is 


Va< z - z > 


■ ///, 


3 

[ l 

i ,J=1 


o i J(z)e i J(z)]dn 


[ 2 . 6 ] 


Even though shape design variable, which is the shape of the domain fl, is 
the only design variable in this case, subscript u is left for general 
treatment. The load linear form [2] of the three dimensional elastic 
solid is 


[2.7] 


* = ///q t l < r1 z 1 ]da + // 2 [ I T^jdr 

u *“ “ 1-1 r 1-1 

0 1 2 

where r , r , and r are clamped, traction free, and loaded boundaries, 
respectively, f = [f*, is the body force, and T « [T 1 , T^, T"*]^ 

is the traction force. 

The variational equation of the three dimensional elastic solid is [2] 

a u ^(z,T) s Vatt. for a11 z [2.8] 

where Z is the space of kinematically admissible displacements; i.e., 

Z = {z G [H 1 ( «) ] 3 : (x) = 0, i = 1, 2, 3, x e r 0 } [2.9] 

For plane elasticity problems in which either all components of stress 
in the x 3 -direction are zero or all components of strain in the x 3 - 
direction are zero, Eq. 2.8 remains valid, with limits of summation 
running from 1 to 2 and an appropriate modification of the generalized 
Hooke's Law of Eq. 2.5. 


C. PLATE/PLANE LEASTIC SOLID 

Consider the plate/plane elastic solid component of Fig. 2.4. The 
energy bilinear form [2] of the component is 

a 

VaUJ) = II DU)[(w n + vw 22 ) w ll + ( w 22 + vw ll ) w 22 

2 . . . 

+ 2(l-v)w 19 w 19 ldQ + // t[ l o 1J (v)e 1J (v)]dfi [2.10] 

U 12 a i ,j=l 

1 p T ~ 

where z * [w, v , v ] is the displacement field. In Eq. 2.10, D(t) - 

Et /( 12( 1-v ) ] , v, and t are flexural rigidity, Poisson's ratio, and 

thickness of the component, respectively. Also, o 1J (v) and e 1J (v) are 

stress and strain due to an in-plane displacement field v = [v , v ] , 

respectively. For this component, the conventional design variable is u = 

t(x) and the shape design variable is the shape of the domain fi. The load 

linear form [2] of the component is 


Figure 2.4 Plate/Plane Elastic Solid Component 

2 . . 2 

\j n (T) 55 // + // t I fVjdn + L [l TVjdr [2.11] 

n n i=l r i=l 

where q, f s [f*, and T = [T*, are lateral load, body force, and 

traction force, respectively, as shown in Fig. 2.5. As in the beam/truss 
component, if there are point loads, a Dirac delta measure [2] can be 
used. 



Figure 2.5 External Loads For Plate/Plane Elastic Solid 



The variational equation of the plate/plane elastic solid component is 


[ 2 ] 


a u^(z,z) * for a11 zeZ [2.12] 

9 19 

where Zc H ( n) x [H ( ft) ] and elements of Z satisfy kinematic boundary 
conditions. For plane leastic solid, kinematic boundary condition is 

v 1 (x) =0, i = 1,2, x 6 r° [2.13] 

For plate, the boundary can be clamped, simply supported, or free edge. 
While the calculation may not be as simple as in the case of beam, the 
variational Eq. 2.12 is valid for all boundary conditions considered [2]. 

Note that Eqs. 2.3, 2.8, and 2.12, the variational equations for 
different structural components are all in the same form. 


3. MATERIAL DERIVATIVE FOR SHAPE DESIGN SENSITIVITY ANALYSIS 

The first step in shape design sensitivity analysis is development of 
relationships between a variation in shape of a structural component and 
the resulting variations in functionals that may arise in the shape design 
problems. Since the shape of domain a structural component occupies is 
treated as the design variable, it is convenient to think of a as a 
continuous medium and utilize the material derivative idea of continuum 
mechanics. In this section, the definition of material derivative is 
introduced and several material derivative formulas that will be used in 
later sections are derived. 

Consider a domain a in one, two, or three dimensions, shown 
schematically in Fig. 3.1. Suppose that only one parameter t defines the 
transformation T, as shown in Fig. 3.1. The 
mapping T : x x^(x) , x e fi, is given by 


Figure 3.1 One Parameter Family of Mappings 


The process of deforming a to by the mapping of Eq. 3.1 may be 
viewed as a dynamic process of deforming a continuum, with x playing the 
role of time. At the initial time x * 0, the domain is a. Trajectories of 
points x e ft, beginning at x * 0, can now be followed. The initial point 
moves to = T(x,x). Thinking of x as time, a design velocity can be 
defined as 


V(x T ,x) = 



3T(x,x) 

3x 


[3.2] 


In a neighborhood of x = 0, under reasonable regularity hypotheses 

T(x,x) = T(x,0) + (x,0) + 0(x 2 ) 

■ x + xV(x,0) + 0(x 2 ) 


Ignoring higher order terms, 

T(x,x) * x + xV(x) [3.3] 

where V(x) = V(x,0). In this paper, only the transformation T of Eq. 3.2 
will be considered, the geometry of which is shown in Fig. 3.2. 

Variations of the domain ft by the design velocity field V(x) are denoted 
as * T(ft,x) and the boundary of ft^ is denoted as r . Henceforth in the 
paper, the term “design velocity" will be referred to simply as 
"velocity". 

1 / 

Let ft be a C -regular open set; i.e., its boundary r is closed and 
bounded and can be locally represented by a C k -function. Let V(x) E R n 
in Eq. 3.2 be a vector defined on a neighborhood U of the closure U of ft 




Figure 3.2 Variation of Domain 


and V(x) and its derivatives up to order k > 1 be continuous. With these 
hypotheses, it has been shown [20] that for small t, T(x,t) is a 
homeomorphism (a one-to-one, continuous map with a continuous inverse) 
from U to U T = T(U,t) and that T(x,t) and its inverse mapping T~*(x x ,x) 
have C k -regularity and has C k -regularity. 

Suppose z x (x T ) is a smooth solution of the elasiticity equations. 

Then the mapping z t (x t ) = z x (x + tV(x)) is defined on 0 and z x (x x ) depends 
on t in two ways. First, it is the solution of the boundary-value problem 
on Second, it is evaluated at a point x^ that moves with t. The 
pointwise material derivative (which is shown to exist in Ref. 2) at 
x e n is defined as 


z(x) 


37 * t (x+tV(x) )| t=0 * 


1 im 

T+0 


z T (x + xV(x) ) - z(x) 

T 


[3.4] 


If z^ has a regular extension to a neighborhood U x of « x , then 

z(x) = z'(x) + Vz T V(x) [3.5] 


where 


z'(x) 


lim 

T+0 


Z x (x) - z(x) 

T 


[3.6] 


is the partial derivative of z. 


One attractive feature of the partial derivative is that, with 
reasonable smoothness assumptions, it commutes with the derivatives with 
respect to [2]; i .e. , 


3z 1 1 
3x ' 


i 


3Xi 


(z ) 


1,2,3 


[3.7] 


A pair of technical material derivative formulas that are used 
throughout the remainder of the paper are summarized in this section. 
Their proofs are presented in Ref. 2. 


Lemma 3.1 : Let i|»^ be a domain functional, defined as an integral over 


♦ l * //„ %<*,> C3.S3 

T 

where f is a regular function defined on ft . If ft has (^-regularity, 

T T 

then the material derivative of ^ at ft is 

= JJ a f * (x)dft + J r f(x) (V T n) dr [3.9] 

or, equivalently, 

i '[ u fj a [f 1 (x) + Vf (x) T V(x) + f(x)div V(x) ]dft [3.10] 

It is interesting and important to note that only the normal component 
(V^n) of the boundary velocity appearing in Eq. 3.9 is needed to account 
for the effect of domain variation. In fact, it is shown by Theorem 3.5.2 
of Ref. 2 that if a general domain functional «|/ has a gradient at ft and 
if ft has C^-regularity, then only the normal component (V^n) of the 
velocity field on the boundary is needed for derivative calculations. 

In contrast to Eq. 3.9, use of the mathematically equivalent result 
given in Eq. 3.10 requires that the velocity field V(x) be defined 
throughout the domain ft. Of course, it must be consistent with 
(V T n) on r. Nevertheless, there are an infinite number of velocity fields 
that satisfy this condition, for each of which the result of Eqs. 3.9 and 
3.10 must be the same. 

Next, consider a functional defined as an integration over r , 


[3.11] 


*2 * It 9 x (x t> dr T 

T 

Lemma 3.2 ; Suppose g T in Eq. 3.11 is a regular function defined 
on r T . If 51 is C k+1 regular, the material derivative of ^ is 

^2 = /p ts'U) + ( v S Tn + H 9( x ) ) (V T n) ] dr [3.12] 

where H is the curvatue of r in and twice the mean curvature in R^. 

4. ADJOINT VARIABLE FORMULATION OF SHAPE DESIGN SENSITIVITY ANALYSIS 

As seen in Section 3, the static response of a structure depends on 
the shape of the domain. Existence of the material derivative z, which is 
proved in Ref. 2, and material derivative formulas presented in Section 3 
are used in this section to derive an adjoint variable method for design 
sensitivity analysis of several functionals. Since the finite element 
method is used for numerical analysis of the structural systems in this 
paper, only the domain method of shape design sensitivity analysis is 
presented in this section. 

The variational equations of several structural components of Eqs. 

2.3, 2.8, and 2.12 on a deformed domain, is of the form 

Vfl < 2 t> *t> 5 fitter- z t> da T 

= llo fTz T da t + VVx s Va for 3,1 z t e Z T 

t r * x 

T 

[4.1] 

where is the space of kinematically admissible displacements on 0^ and 
c(*,*) is a bilinear mapping that is defined by the integrand of Eqs. 2.1, 
2.6, and 2.10. 

Taking the material derivative of both sides of Eq. 4.1, using Eqs. 
3.10 and 3.11 and noting that the partial derivatives with respect 
to r and x commute, 

^ a u,^ 2 ’ 2 ^ 5 a u,V^ z ’^ + a u,Q (z ’ z) * \j,V^’ f0r 311 ~ ze 1 t4 * 2 - i 


where, using Eq. 3.5, 


[ a u , n ( z »z)] = // n [ c ( z * z ) + c(z',z) + Vc(z,z) V + c(z,z)div V] da 
= //Jc( 2 ,z - Vz T V) + c(z - Vz T V,z) + Vc(z,7) T V 
+ c(z,7)div V] da [4.3] 


and 


i u , v (i) = //q [f T i + v(f T z) T V + f T zdiv V]da 

+ J *> [T T z'+ (v(T T z) T n + HT T z)(V T n)]dr 
r 

= ff Q [f T (z - Vz T V) + v(f T z) T V + f T zdi v V]da 

+ I p [T T (z - Vz T V) + (v(T T z) T n + HT T z)(V T n)]dr [4.4] 

r 

The fact that the partial derivatives of the coefficients, which depend on 

cross-sectional area and thickness, in the bilinear mapping c(*,») are 

zero has been used in Eq. 4.3 and f' = T‘ = 0 has been used in Eq. 4.4. 

0 12 

For boundary variations, it is supposed that the boundary r = r U r U r 

is varied, except that the curve 31^ that bounds the loaded surface is 

2 

fixed for three dimensional elastic solid, so the velocity field V at 3r 

? 

is zero. For the case in which 3r is not fixed, variation of the 

2 

traction term in Eq. 4.1 (given as an integral over r ) gives an 
additional term that was not discussed in lemmas. For this case, the 
interested reader is referred to Ref. 21. For plane elastic solid 
component case, these additional terms will be given in Section 5. 

For z x , select z t (x + tV(x)) = z(x); i.e., choose z as constant on the 
line x t = x + tV(x). Then, since H m (a) is preserved by T(x,t) 
(homeomorphism property noted in Section 3), if 5 is an arbitrary element 
of H m (a) that satisfies kinematic boundary conditions on r, z t is an 
arbitrary element of H m (a T ) that satisfies kinematic boundary conditions 
on r . In this case, using Eq. 3.5, 

1 - 1 -T 

z = z +VzV = 0 


[4.5] 



From Eqs. 4.2, 4.3, and 4.4, using Eq. 4.5, 


a u,V (z,i) = //n [-cU»’z T V) - c( Vz T V,z) 
+ vc(z,z) T V + c(z,z)div Vjdn 


C4.6] 


and 


* u , v (2) = // n [z T (vf T V) + f T zdiv V]dn 

+ / ? [-T T (Vz T V) + (v(T T z ) T n + HT T z)(V T n)]dr [4.7] 

r 

Then, Eq. 4.2 can be rewritten to provide the result 

a u>n (z»z) 55 * u>v (z) ' a u,V^ z,z ^ for a11 zeZ [4.8] 

Consider a displacement functional that defines the displacement at 
nodal point x e ft 

i|>. = z(x) = // 6(x-x)z(x)dn [4.9] 

* 

where S(x) is the Dirac delta measure at the origin. Taking the first 
variation of Eq. 4.9, using the material derivative, 


= z(x) = // n 6(x-x)z(x)dfl [4.10] 

The objective now is to obtain an explicit expression for ^ in terms 
of the velocity field V, which requires eliminating z. An adjoint 
equation is introduced by replacing z E Z in Eq. 4.10 by a virtual 
displacement IeZ and equating terms involving X to the energy bilinear 
form, yielding the adjoint equation for the adjoint variable X, 

a u fi (X,X) = Jf n 6(x-x)X(x)dO, for all XGZ [4.11] 

Denote the solution of Eq. 4.11 as X^. 

To take advantage of the adjoint equation, evaluate Eq. 4.11 at 
X = z, since z e Z [2], to obtain the expression 


[4.12] 


a u,Q (x{1) * i) = ^ft ^ x ' x )^ x ) dn 


(1) 


Similarly, evaluate the identity of Eq. 5.8 at z = X' , since both are in 
Z, to obtain 




[4.13] 


Recalling that the energy bilinear form a -(•»•) is symmetric in its 

U ))■ 

arguments, the left sides of Eqs. 4.12 and 4.13 are equal, so 


// n 6(x-x)z(x)dft = A^ V (X^^) - a u >v (z,^^) 


[4.14] 


Using Eqs. 4.14, Eq. 4.10 yields 




[4.15] 


Explicit expressions of the terms in Eq. 4.15, for each structural 
component can be obtained using Eqs. 4.6 and 4.7. These explicit 
expressions will be derived in Section 5. This order of presentation was 
chosen to show basic idea of the adjoint variable method without 
complicate derivation of expressions. 

Note that evaluation of the design sensitivity formula of Eq. 4.15 
requires solution of Eq. 4.1 for z. Similarly, Eq. 4.11 must be solved 
for the adjoint variable A^. This is an efficient calculation, using 
finite element analysis, if the boundary-value problem for z has already 
been solved, requiring only evaluation of the solution of the same set of 
finite element equations with different right side, called an adjoint 
load. 

Next, consider a locally averaged stress functional over a test 

volume ft of the three dimensional elastic solid, 

P 

/// Q g(o(z))dn 

*2 = Ma 9(°( z ) ) m p dft " fj '~~ da *- 4,16 - 1 

P 

where o denotes the stress tensor, ftp is an open set, and mp is a 
characteristic function that is constant on ft p , zero outside of ft p , and 



whose integral is 1. Here, g is assumed to be continuously differentiable 
with respect to its arguments. Note that g(a(z)) might involve principal 
stresses, von Mises failure criterion, or some other material failure 
criteria. Taking the first variation of Eq. 4.16, using Eq. 3.10 [10], 


*2 = t /// 0 (9* + v 9 Ty + 9 dl ' v V)dn/// fl d 0 
P P 

- JJJ n 9 da//J fi div Vd«]/ C/// n dn) 2 

P P P 

= /// o I 9 i i (z) [o 1J (z) - a 1J (Vz T V)]m da 

n i ,j=l o ,J p 

+ ffl a nig i i ( 2 )°k j ( z ) vk ]m D dO + ffl a 9 div Vm da 
“ k=l i ,j=l a J k p w p 

- J/J fi 9 m p dJ2 ///fl m p div Vda [4.17] 


It can be shown that 

a ij (Vz T V) = l C^ k£ (Vz^ V + Vz k VJ 
k ,£=1 1 K 


and 


[4.18] 


3 3 . j 

I oj j (z)V k = I C ijla (Vz k V) [4.19] 

k=l K k,i»l 

Using these results, Eq. 4.17 becomes 
*2 “ J//J I 9 ij (z)o lj (z)]m do 

-/// Q l [ l 9 ii (z)C ijk4 (Vz kT V )]m da 
“ i,j=l k,*=l a 1J * P 

+ fff a 9 div Vm p da - J// Q gm p da fff Q m p div Vda [4.20] 

As in the displacement functional case, an adjoint equation is 
introduced by replacing z eZ in the term on the right of Eq. 4.20 by a 
virtual displacement ie Z and equate the result to the energy bilinear 
form. 



for all X e Z 


[4.21] 


* n (X ’*) s // t l 9 ij (z)o 1 J (X)]m da, 

u,n Q i ,j=l a J p 


Denote the solution of Eq. 4.21 as X^* By the same method used for the 
displacement functional, the sensitivity formula is obtained as 


♦j - u <* (2) >. - < 


( 2 ) 


u.V 


u ,V 


-III It! 9 11 (z)C 1Jk *(vi kT V )]m d# 

a i ,j=l k ,t=l o 1J * p 

+ /// 9 div Vm da - /// gm da /// m div Vda 


[4.22] 


a 


a 


a 


where explicit expressions of the first two terms in Eq. 4.22 for the 
three dimensional elastic solid can be obtained using Eq. 2.6, 2.7, 4.6, 
and 4.7. These explicit expressions will be derived in Section 5. Note 
that these terms have the same form as those of Eq. 4.15 for the three 
dimensional elastic solid. The difference is that terms in Eq. 4.15 are 
evaluated at X^ and terms in Eq. 4.22 are evaluated at X^. That is, 
once the expressions for terms in Eq. 4.15 are derived, they can be used 
for different functionals. 

Finally consider a locally averaged stress functional over a test 
area in a plate/plane elastic solid component. 


h = // [ 9 1 (t,w ij ) + g 2 (a(v))]m n da 


[4.23] 


1 2 

where g i (t,w^j) and g (a(v)) are principal stress, von Mises yield stress, 
or some other stress measures due to lateral displacement w and in-plane 
displacement field v, respectively. Here, g^(t,w 1 -j) is measured at the 
extreme fiber and nip is a characteristic function on that is constant 
on a^, zero outside Op, and whose integral is 1. Taking the first 
variation of Eq. 4.23, using Eq. 3.10, 


*3- fi t“ij - < 5 '" Tv >i J ] v a 


a i ,j=l i j 


+ // [ I g 2 ij - o lj (v)]m n da + J/ di v(g 1 V) m n da 


a i,j=l a 


a 


[4.24] 



2 


l 

U = 1 


[ l g 2 , i (v)C 1J ’ k4 (Vv k V )] id q 
k,£»l a J * p 


+ // g 2 div V m da - // ( g^+g 2 ) m dil // m div V da 

a p a p a p 

Define an adjoint equation by replacing w and v in Eq. 4.24 by virtual 
displacements "n and T, respectively, and equate terms involving n and £ in 
Eq. 4.24 to the energy bilinear form, 

2 , 

a u 8 //[ l 9 W n- Jm dn 

u ’“ o i,j=l w ij 1J p 

2 2 ii_ 

+ //[ l 9 -h a J (C)]mdn, for all XGZ [4.25] 

0 i ,j=l o 1J p 

1 2 T 

where X = [n,5 ,6 ] is an adjoint variable. Denote the solution of Eq. 
4.25 as x^. By the same method used for the displacement functional, 
the sensitivity formula is obtained as 

2 

*3 * *u v (x(3)) " K v (z ’ x(3)) " // ( l gj (Vw T V)]m dn 
J u.v u,v oi,j=l w ij p 

2 2 . . T 

+ //div(g 1 V)m dfl - // l [ l g 2 1i (v)C ijkA (W k V )]m do 
a p ft i,j=l M=1 o 1J * p 

+ // g 2 div V m dfl - // (g*+g 2 )m do // m div V dfl [4.26] 

a p ft p fl p 

where explicit expressions of the first two terms in Eq. 4.26 for the 
plate/plane elastic solid component can be obtained using Eqs. 2.10, 2.11, 
4.6 and 4.7. These explicit expressions will be derived in Section 5. As 
in the displacement functional case, evaluation of the design sensitivity 
formula of Eq. 4.26 requires solutions z and X^ of Eqs. 4.1 and 4.25. 
Design sensitivity information for locally averaged stress functional over 
a test length Op in a beam/truss component can be derived using the same 
procedure. 


5. SHAPE DESIGN SENSITIVITY ANALYSIS OF STRUCTURAL COMPONENTS 


In this section, explicit expressions for terms in Eqs. 4.6 and 4.7 
are derived for each structural components by identifying bilinear 
mapping c(*,*) and loading terms. The result is standard expressions that 
can be used for design sensitivity analysis of different functionals. 

These results can also be used for design sensitivity analysis of built-up 
structures which will be shown in Section 9. 


A. BEAM/TRUSS 

Using energy bilinear and load linear forms of Eqs. 2.1 and 2.2 for 
beam/truss component, Eqs. 4.6 and 4.7 become 

£ 

a 1 MJ) * / {-EI*[3w* w* V + (w 1 w 1 + )V ] 
u,V v ' 1 1 xx xx x v xx x x xx' xx J 

£ 

+ EI x«x V > d * ♦ / l- El2 [ 3 *4"xx V x 
+ (w 2 ^ + w 2 ^ 2 )V ] + EI 2 w 2 ^ V} dx 

' XX X X XX / XX J X XX XX J 

+ n- GJ w x + GJ xW> dx 


* l (- hEv x v x V x + h x Ev x v x V) dx 


[ 5 . 1 ] 


and 


t; #v (z) = / (q^V + q 1 w 1 V x ) dx + / (q^V + q 2 w 2 V x ) dx 

l _ _ £ _ 

+ / (r 0V + rev ) dx + / (f Y vV + fvV ) dx 


[5.2] 


B. THREE DIMENSIONAL ELASTIC SOLID 

Using energy bilinear and load linear forms of Eqs. 2.6 and 2.7 for 
three dimensional elastic solid, Eqs. 4.6 and 4.7 become 


and 


a u V (z,z) = " Mo X [<J lj (z)e 1J (Vz T V) + a 1J (z)e 1J (Vz T V)]dQ 
U,V " i ,J = 1 

+ /J7fl v [ I a 1J (z)e lj (7)] T Vdfl + /// [ £ cr 1 J (z) e 1 J (7) ]di v Vdfl 

i,j=l i,j=l 

[5.3] 


3 -i i T 

/ nfl 


3 i-i i 


A u V (z) “ ///ft I z < vf v ) dn + /// Q [ I f 2 ] div Vdn 

u,v i=1 “ i=l 


3 .. ,T 3 . , T 3 . , 


+ // 2 {- I T^z 1 V) + (V[ l T^z 1 ] T n + H[ £ T'z 1 ] )(V 'n)}dr 

r i=l i=l i=l 


1=1n/»T- 


[5.4] 


It can be verified that 

l a lj (z)e 1J '(v7 T V) = l a lj (z)(v7 1 . T V + v7 lT V.) 
l.j-1 U=1 0 J 

and 


[5.5] 


V[ l a ij '(z)e ij (7)] T V * l [a ij (z)(v7fv) + a ij (7) ( Vzf V) ] [5.6] 

l.j-1 l.j-1 J J 

1 2 3 T 

where V. = [V-,V.,V.] . Using the above results, Eqs. 5.3 becomes 

J J J J 

3 .. .T .. .T 

J / 7 U M \ 4- rr^Jf-yWUT^ \J ) j(jft 


^ jV (z»z) = " ///ft . I i a J ( z )( Vz v j) + a J U)( Vz Vj.; 

i >j = l 

+ J/JJ I o 1J (z)e 1J (7)] div Vdfi 
i,j=l 


[5.7] 


C. PLATE/PLANE ELASTIC SOLID 


As in the beam/truss and three dimensional elastic solid components, 
explicit expressions for terms in Eqs. 4.6 and 4.7 can be obtained using 
energy bilinear and load linear forms of Eqs. 2.10 and 2.11 for 
plate/plane elastic solid component. For plane elastic solid component, 
Eqs. 5.7 and 5.4 remain valid, with limits of summation running from 1 to 
2 and an appropriate modification of the generalized Hooke's Law of Eq. 
2.5. Equations 4.6 and 4.7 become, for plate/plane elastic solid 
component. 


a u,V^ Z * Z ^ ~ ^ ( w ii w n^i + w 22 w 22^2^ + t v ( w u w 22 + W 22 W 1 1 ^ 

' ( w ll"ll + * 22 * 22 ^ + 2 ( 1 - V ) w 12 w 12 ] div v 
+ 2 Kl"l2 + W 12*ll + *22*12 + W 12*22^ V 2 + V l^ 

+ [ Wl w u + + v (w x w 22 + w 22 w 1 )]v} 1 

+ [w^w 22 + * 22*1 v ^*l*ll + *11*1^^22 

+ t w 2"ll * w n"2 * u ("2“22 + 'V^ll 

+ [w 2 w 22 + *22*2 + v ( w ? w i i + w i i w t? 3 ] ^ 


2 11 11 2 ^ 22 
77 Vfl . r . . — ... 77 Vi2 


+ 2(l-v) [(wjWj 2 + W 12 W 1 ^12 + ( W 2 W 12 + *12*2^12^ 


Et 


[wj i W-| i + v ( w n w 22 + W 22 W 11 ^ + w 22 w 22 


4(l-v 2 ) 1 11 11 

+ 2(l-v)Wj 2 Wj 2 ]vt^V dfl 

2 .. .T .. .T 

+ II {-t l [o 1J (v) (Vv 1 V.) + o 1J (v) (Vv 1 V.) 
n i ,j=l J J 

2 

- (v) (V)divV ] + l o^vJ.^ttMda 

i ,j»l 


[5.8] 


and 


*[, >v (z) = H [w(Vq T V) + qw div V]dfl 


+ // l [7 1 (Vf 1 V) + fV div V]d n 
fl i = l 

+ I l {-T 1 (w 1 V) + [v(tV) n + HT 1 v 1 ] (V T n ) }dT 
r i=l 


+ 


2 


l 

1=1 


[tVv 


Tp- 


- tVv 


T Pi 


[5.9] 


2 

In Eq. 5.9, H is the curvature of the loaded boundary r and the last two 
terms on the right account for corner effects due to movement of points pi 



and P2 in Fig. 2.5 [21]. In these two terms, the notation p.. indicates 
that the terms are evaluated at point p^ and Vj is the component of 
velocity V tangent to r, which is positive if it is in a counter-clockwise 
direction in Fig. 2.5. 

Results given in Eqs. 5.1, 5.2, 5.4, and 5.7 - 5.9 can be used in Eqs. 
4.15, 4.22, and 4.26 for each structural component and each functional. 
This allows one to develop a modular computer program that will carry out 
numerical integrations of terms in Eqs. 5.1, 5.2, 5.4, and 5.7 - 5.9 using 
the same shape functions that are employed in finite element analysis 
codes. The result will then be a general algorithm and numerical method 
for design sensitivity analysis that can be implemented with existing 
finite element codes which will be discussed in Section 6. 

6. IMPLEMENTATION OF DESIGN SENSITIVITY ANALYSIS WITH 

EXISTING FINITE ELEMENT CODES 

To obtain design sensitivity information, Eqs. 4.1 and 4.11 must be 
solved for displacement functionals and Eqs. 4.1, 4.21 and 4.25 must be 
solved for stress functionals. Once the original and adjoint structures 
are solved, one can integrate Eqs. 4.15, 4.22, and 4.26 numerically to 
obtain the desired sensitivity information. The finite element method can 
be viewed as an application of the Galerkin method to Eqs. 4.1, 4.11, 

4.21, and 4.25 for an approximate solution of the boundary-value 
problem. Note that the energy bilinear forms for Eqs. 4.1, 4.11, 4.21, 
and 4.25 are same. Hence, the adjoint structures of Eqs. 4.11, 4.21, and 
4.25 are the same as that of Eq. 4.1, with different adjoint loads. The 
adjoint load of Eq. 4.11, is a simple unit load at the point x in the 

A 

positive direction of z(x). To calculate the adjoint load using the load 
functional on the right side of Eqs. 4.21 and 4.25, one should use the 
same shape functions that are used in the finite element code. Since m p 
is a characteristic function defined on finite element n p , numerical 
integration of the load functional is done on ft p only and the adjoint 
equivalent nodal force acts only on the nodal points of ft p . 

For numerical implementation with existing finite element codes, one 
can proceed as in the flow chart of Fig. 6.1. In the beginning, the model 
is defined by identifying the finite element model, original structural 
load, design variables, and constraint functionals. In the next step, an 




Figure 6.1 Flowchart of Design Sensitivity Calculation Procedure 












existing finite element code is called to obtain structural response. 

With the structural response obtained, one calculates an adjoint load, 
external to the finite element code, using the shape functions of the 
code. The adjoint load is then input to the finite element code, to 
obtain an adjoint response for each constraint functional. For adjoint 
analysis, one can use the multiloading (restart) option of the finite 
element code, so that only forward and backward substitutions are 
performed to obtain each adjoint response. Using the original and adjoint 
structural responses, design sensitivity information is calculated for 
each constraint functional, by carrying out only numerical integration. 
This procedure allows one to carry out calculations outside finite element 
codes, using postprocessing data only. That is, the design sensitivity 
software does not have to be imbedded into finite element codes. 

Moreover, the method does not require differentiation of stiffness and 
mass matrices and the uncertainty of numerical accuracy associated with 
selection of a finite difference perturbation can be eliminated. 


7. NUMERICAL EXAMPLES 

Substantial numerical experimentation has been carried out using the 
material derivative shape design sensitivity analysis formulation, with 
the boundary method. Good results have been reported [2,23] for a variety 
of single structural components. These studies have shown that great care 
must be taken in projecting stress information to the boundary to achieve 
acceptable design sensitivity accuracy. Higher order elements and 
extrapolation from Gauss points have been shown to be essential in 
achieving acceptable accuracy. Substantially inaccurate results have been 
observed when low order elements are used and elementary boundary 
projection approaches are employed. 

Numerical experimentation with the domain method [15,16,18,22] has 
indicated consistently good results for structural components, without the 
requirement for sophisticated elements, clever boundary projection 
methods, or drastically refined grids. In order to be more quantitative, 
two examples are discussed to permit numerical comparison. 

Consider a plane elastic solid that is composed of two materials of 
substantially different modulus of elasticity (E /E = 7.65) and subjected 
to simple tension, as shown in Fig. 7.1. The finite element configuration. 




40 N/cm 


Figure 7.1 Interface Problem 

dimensions, material properties of each body, and loading conditions are 
shown in Fig. 7.1. Body i occupies domain fl 1 , i=l,2, AB is the interface 
boundary y, and and are the clamped and loaded boundaries, 
respectively. The design variable b controls the position of the inter- 
face boundary y» while the overall dimensions of the structure are fixed. 

The expression for design sensitivity of the von Mises yield stress 
functional associated with interface boundary movement with the domain 
method is obtained by simply adding results of Eq. 4.26 for both segments 
of the structure. For the plane stress interface problem, terms in Eq. 
4.26 due to plate bending must be dropped. For the boundary method, 
design sensitivity computations are carried out in Ref. 10 (Eq. 42) that 
is analytically equivalent to the result of the domain method. 

For numerical computation, the finite element method is used to 
approximate the state and adjoint equations of Eqs. 4.1 and 4.25, 
respectively. In order to compute the design sensitivity expressions of 
Eq. 4.26 one must define a design velocity field V that satisfies 
regularity properties defined in Refs. 2 and 9, in terms of variations in 
the design variable b. To have a continuous design velocity field, one 
may define 



[7.1] 



The finite element model shown in Fig. 7.1 contains 32 elements, 121 
nodal points, and 224 degrees-of-freedom. The 8-noded isoparametric 
element is employed for design sensitivity analysis. For the boundary 
method, stresses and strains are obtained at Gauss points and extrapolated 
to the boundary to obtain accurate results on the boundary [23]. Define 

i 2 

t and i|> as the functional values for the initial design b and modified 

design b + 6b, respectively. Let Ai|> = i|> - ^ and let be the predicted 

difference from sensitivity analysis. The ratio times 100 is used 

as a measure of accuracy; i.e., 100% means that the predicted change is 

exactly the same as the actual change. Notice that this accuracy measure 

will not give meaningful information when Aiji is very small compared 

to ^ , because the difference a<j> may lose precision due to the 
2 1 

subtraction \|> - . 

Numerical results with a 3% design change; i.e., 6b = 0.03b, are shown 
in Table 7.1 for the boundary method and in Table 7.2 for the domain 
method. Due to symmetry, sensitivity results for only the lower half of 
the structure are given. These results indicate that the domain method 
gives excellent results, whereas accuracy of the boundary method is not 
acceptable. For elements 22 and 29, the predicted values are less 
accurate than others. However, the magnitude of actual differences A<|> for 
those elements are smaller than others, so Aij> may lose precision. 

A disadvantage of the domain method is that a velocity field must be 
defined in the domain and satisfy regularity properties. There is no 
unique way of defining domain velocity fields for a given normal velocity 
field (V^n) on the boundary. Also, numerical evaluation of the 
sensitivity result of Eq. 4.26 is more complicated than evaluation of Eq. 
42 of Ref. 10, since Eq. 4.26 requires integration over the entire domain, 
whereas Eq. 42 of Ref. 10 requires integration over only the variable 
boundary. This problem can be alleviated by introducing a boundary layer 
[16] of finite elements that vary during the perturbation of the shape of 



Table 7.1. Boundary Method for Interface Problem 


El. 

No. 


♦ 2 

Aip 


(*'/A*xlO0)X 

1 

393.01304 

393.17922 

0.16618 

0.20403 

122.8 

2 

364.37867 

364.76664 

0.38796 

0.67218 

173.3 

5 

388.07514 

388.36215 

0.28701 

0.56684 

197.5 

6 

402.26903 

402.83406 

0.56503 

0.42080 

74.5 

9 

386.43461 

386.84976 

0.41515 

-0.08520 

-20.5 

10 

407.14612 

407.48249 

0.33637 

0.14159 

42.1 

13 

388.59634 

388.95414 

0.35780 

-0.53089 

-148.4 

14 

379.04276 

379.25247 

0.20971 

-1.90134 

-906.6 

17 

441.68524 

442.25032 

0.56507 

-13.85905 

-2452.6 

18 

424.05820 

425.22910 

1.17089 

-13.63066 

-1164.1 

21 

424.19015 

424.70840 

0.51825 

-0.21408 

-41.3 

22 

378.85433 

378.97497 

0.12064 

0.76770 

636.4 

25 

407.71528 

408.23368 

0.51840 

0.49780 

96.2 

26 

387.87307 

387.32342 

-0.54962 

-0.48837 

88.9 

29 

400.61014 

400.60112 

-0.00903 

0.01423 

-157.7 

30 

394.61705 

394.00702 

-0.61003 

-0.57794 

94.7 



Table ; 

7.2. Domain Method for Interface Problem 

El. 

No. 


* 2 

Alp 

V 

(f /A*xl00)% 

1 

393.01304 

393.17922 

0.16618 

0.17954 

108.0 

2 

364.37867 

364.76664 

0.38796 

0.37840 

97.5 

5 

388.07514 

388.36215 

0.28701 

0.28671 

99.9 

6 

402.26903 

402.83406 

0.56503 

0.59634 

105.5 

9 

386.43461 

386.84976 

0.41515 

0.41515 

100.6 

10 

407.14612 

407.48249 

0.33637 

0.33637 

109.6 

13 

388.59634 

388.95414 

0.35780 

0.37548 

104.9 

14 

379.04276 

379.25247 

0.20971 

0.20159 

96.1 

17 

441.68524 

442.25032 

0.56507 

0.57069 

101.0 

18 

424.05820 

425.22910 

1.17089 

1.12871 

96.4 

21 

424.19015 

424.70840 

0.51825 

0.53919 

104.0 

22 

378.85433 

378.97497 

0.12064 

0.06396 

53.0 

25 

407.71528 

408.23368 

0.51840 

0.51710 

99.7 

26 

387.87307 

387.32342 

-0.54962 

-0.56083 

102.0 

29 

400.61014 

400.60112 

-0.00903 

-0.00298 

33.0 

30 

394.61705 

394.00702 

-0.61003 

-0.58529 

95.9 


a structural component. This approach is illustrated schematically in 
Fig. 7.2. The domain ft is divided into subdomains 8^ and with inner 
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Figure 7.2 Boundary Layer 

core held fixed and only boundary layer modified. In this way, the 
velocity field need be defined only on ^ The thickness of the boundary 
layer ^ will depend on trade-offs between numerical accuracy and 
numerical efficiency. In practice, ^ can be a substructure of the finite 
element model. 

To demonstrate feasibility of the boundary layer approach, two 
examples are solved by the boundary layer approach. The first example is 
the plane stress interface problem discussed in this section. For a body 
of given geometry there is a large number of possible boundary-layers, 
some of which are better than others, from the viewpoint of accuracy and 
efficiency. It is difficult to estimate the size and location of the best 
boundary-layers in advance. They can be determined by analyzing the 
structure and measuring the strain energy density [24]. 

The boundary-layer is chosen to include elements 13 thru 20 in Fig. 
7.1. The design variable b for this case is distance between node 51 and 
node 65 in Fig. 7.1. Consequently, regions outside the boundary-layer 
remain unchanged. Numerical results with a 3% design change are shown in 
Table 7.3 for the boundary-layer approach. Due to symmetry, the shape 
design sensitivity analysis results of the lower half of the structure are 
shown. Shape design sensitivity analysis results obtained with the 
boundary-layer approach are excellent, as shown in Table 7.3. 



Table 7.3. Boundary Layer Approach for Interface Problem 
(E2/E 1 = 7.65) 


El. 

No. 

+ 1 

* 2 

Aip 

V 

( ip* /Atpxl00)% 

1 

393.01304 

393.07967 

0.06663 

0.06770 

101.6 

2 

364.37867 

364.29542 

-0.08325 

-0.08412 

101.0 

5 

388.07514 

388.20344 

0.12830 

0.12916 

100.7 

6 

402.26903 

402.19633 

0.07270 

-0.07126 

98.0 

9 

386.43461 

386.47687 

0.04227 

0.04103 

97.1 

10 

407.14612 

407.46571 

0.31960 

0.32954 

103.1 

13 

388.59634 

388.70300 

0.10666 

0.10834 

101.6 

14 

379.04276 

379.52487 

0.48210 

0.48590 

100.8 

17 

441.68524 

442.17008 

0.48484 

0.47615 

98.2 

18 

424.05820 

425.31717 

1.25897 

1.23636 

98.2 

21 

424.19015 

424.44270 

0.25254 

0.25680 

101.7 

22 

378.85433 

378.99459 

0.14025 

0.11902 

84.9 

25 

407.71528 

407.99105 

0.27577 

0.27248 

98.8 

26 

387.87304 

387.59778 

-0.27526 

-0.27540 

100.0 

29 

400.61014 

400.62571 

0.01557 

0.01543 

99.1 

30 

394.61705 

394.45461 

-0.16244 

-0.16104 

99.1 

Next, to test validity of the 

boundary-layer approach 

, Young's modulu 

is changed to E 1 = 0 

.2 MPa and E 2 

= 100 MPa for fl* and S? 

, respectively. 

In other words, the 

ratio between 

E 2 and E* is raised to 

500, from 7.65, 


to check a more severe condition. Design sensitivity results are given in 
Table 7.4. Accuracy of design sensitivity is excellent. For elements 9 
and 22, the magnitude of actual change are small, so finite differences 
may not be accurate. Numerical results obtained with the boundary method 
given in Ref. 16, indicates that worse results arise if the ratio E 2 /E 1 is 
increased. 

Results for the plane stress interface problem clearly indicate that 
the boundary approach may have considerable difficulty in handling 
problems with singular characteristics. Accuracy of the boundary approach 
rapidly deteriorates in the vicinity of a singularity. On the other hand, 
the boundary-layer approach can give good sensitivity results throughout 
the domain. Also, in this interface problem, 56% of cpu time is saved by 
using the boundary-1 ayer approach instead of the domain approach, without 
sacrificing accuracy of design sensitivity. 

Next, the classical fillet shown in Fig. 7.3 is used to study accuracy 
of the boundary layer approach. The design for this problem is the shape 


Table 7.4. Boundary Layer Approach for Interface Problem 
(E 2 /E 1 = 500) 


El. 

No. 

t 1 

'P 2 

A<|> 

V 

(*'/A*xl00)X 

1 

392.30557 

392.39359 

0.08802 

0.08926 

101.4 

2 

365.24352 

365.14384 

-0.09967 

-0.10047 

100.8 

5 

386.63113 

386.77199 

0.14086 

0.14138 

100.4 

6 

402.93637 

402.89886 

-0.03752 

-0.03477 

92.7 

9 

386.58037 

386.56330 

-0.01706 

-0.02016 

118.2 

10 

403.05338 

403.58214 

0.52876 

0.54150 

102.4 

13 

392.16507 

392.15462 

-0.01046 

-0.01033 

98.8 

14 

365.55409 

366.17638 

0.62229 

0.62552 

100.5 

17 

477.04122 

478.18542 

1.14420 

1.11943 

97.8 

18 

440.11698 

442.56025 

2.44327 

2.39258 

97.9 

21 

442.63898 

443.09215 

0.45317 

0.46452 

102.5 

22 

362.56664 

362.67559 

0.10895 

0.07117 

65.3 

25 

412.83466 

413.32142 

0.48676 

0.48000 

98.6 

26 

379.90505 

379.41355 

-0.49149 

-0.49033 

99.8 

29 

401.08818 

401.11877 

0.03059 

0.03031 

99.1 

30 

391.19616 

390.92213 

-0.27403 

-0.27158 

99.1 



IOO lb/ in 


of the varying boundary between points A and B, without moving these 
two points. B-spline representation is used for the varying boundary r^. 
Due to symmetry, only the upper half of a fillet is analyzed. Dimensions 
of the structure and applied loads are given in Fig. 7.3. For material 
property. Young's modulus and Poisson's ratio are 3.0 x 10 7 psi and v = 
0.293, respectively. The segment is the center-line of the fillet 





and r 2 is the uniformly loaded edge. Sensitivity of von Mises stress 
averaged over individual finite elements is employed to test accuracy of 
the boundary-layer approach. The expression for design sensitivity is 
obtained from Eq. 4.26 with nonzero velocity field on the boundary layer. 

The boundary-layer (27% of the total area) shown in Fig. 7.4 is chosen 
after analyzing the structure and measuring the strain energy density. In 
Fig. 7.5, a finite element model with optimized boundary profile and 
319 elements and 1994 active degrees-of-freedom is shown. The element 
type used is an 8-node isoparametric element. 


A 



Figure 7.4 Boundary-Layer of Fillet 
A 



Figure 7.5 Finite Element Mesh of Fillet 

In Table 7.5, shape design sensitivity results for a fillet with 
optimized boundary profile (see Fig. 7.5) are given, obtained with 0.1% 
design perturbation. From Table 7.5, it can be seen that this approach 
can yield excellent shape design sensitivity results. 


Table 7.5. Boundary Layer Approach for Fillet 


El. 

No. 

v 1 

* 2 

Aip 


( <p* / Aipxl00)% 

3 

306.91341 

307.18632 

0.27291 

0.27802 

101.9 

13 

326.43889 

326.65762 

0.21873 

0.22298 

101.9 

23 

386.17664 

386.33965 

0.16301 

0.16693 

102.4 

33 

471.26543 

471.38415 

0.11872 

0.12219 

102.9 

43 

570.91518 

570.95727 

0.04209 

0.04294 

102.0 

53 

669.26309 

669.16035 

-0.10274 

-0.10693 

104.1 

63 

736.39264 

736.12922 

-0.26341 

-0.27187 

103.2 

73 

682.49507 

682.26122 

-0.23385 

-0.23826 

101.9 

83 

761.77304 

761.44391 

-0.32913 

-0.33921 

103.1 

93 

911.47561 

911.14248 

-0.33314 

-0.33869 

101.7 

103 

764.32105 

764.00138 

-0.31967 

-0.32424 

101.4 

113 

884.95173 

884.66270 

-0.28902 

-0.28788 

99.6 

123 

768.06301 

767.74016 

-0.32285 

-0.32715 

101.3 

133 

857.37113 

857.10506 

-0.26607 

-0.26249 

98.7 

143 

999.87828 

999.87875 

0.00047 

0.00045 

95.4 

153 

1009.58379 

1009.52243 

-0.06137 

-0.05071 

82.6 

163 

1000.99060 

1000.92960 

-0.06100 

-0.05863 

96.1 

173 

999.42406 

999.33430 

-0.08976 

-0.08873 

98.8 

183 

1001.15587 

1001.04706 

-0.10881 

-0.10766 

98.9 

193 

958.70044 

958.26124 

-0.43919 

-0.44416 

101.1 

203 

980.40747 

980.03216 

-0.37531 

-0.37607 

100.2 

213 

993.57091 

993.33304 

-0.23787 

-0.23708 

99.7 

223 

1000.47920 

1000.29798 

-0.18122 

-0.17969 

99.2 

233 

762.37599 

762.01522 

-0.36076 

-0.40403 

112.0 

243 

778.19389 

777.74914 

-0.44475 

-0.45837 

103.1 

253 

881.29226 

880.57387 

-0.71838 

-0.69488 

96.7 

263 

1220.22663 

1219.47520 

-0.75143 

-0.72158 

96.0 

273 

835.72273 

835.22544 

-0.49729 

-0.48555 

97.6 

283 

922.89834 

922.29381 

-0.60453 

-0.58846 

97.3 

293 

1033.88104 

1033.07352 

-0.80752 

-0.77133 

95.5 

303 

1093.40867 

1090.74702 

-2.66165 

-2.67483 

100.5 

313 

936.38542 

936.01573 

-0.36969 

-0.37920 

102.6 


8. AUTOMATIC REGRIDDING FOR SHAPE DESIGN 

For numerical implementation of shape design sensitivity analysis, one 
must parameterize the boundary r of the domain ft. For this purpose, one 
may use Bezier curves or surfaces [25]. The next step is to develop a 
general method of defining and computing a velocity field in the domain, 
in terms of the perturbation of the boundary r. Moreover, the velocity 
field must satisfy certain regularity conditions. It is shown in Refs. 2 
and 9 that C^-regular and C^-regular velocity fields are sufficient for 



shape design sensitivity analysis of truss and elastic solid problems and 
beam and plate problems, respectively. However, observing Eqs. 5.1, 5.2, 
5.4, and 5.7 - 5.9, one may relax these regularity conditions. That is, 
for truss and elastic solid problems, the highest order derivative of the 
velocity field that appears in Eqs. 5.1, 5.2, 5.4, and 5.7 - 5.9 is one. 
Thus, one may use a C°-regular velocity field with an integrable first 
derivative. Similarly, one may use a C^-regular velocity field with an 
integrable second derivative for beam and plate problems. Therefore, 
regularity of the velocity field must be at least at the level of 
regularity of the displacement field of the structural component 
considered. This suggests use of displacement shape functions to 
systematically define the velocity field in the domain. Moreover, one 
can select a velocity field that obeys the governing equation of the 
structure. That is, the perturbation of the boundary can be considered as 
a displacement at the boundary. With no additional external forces and a 
given displacement at the boundary, one can use the finite element code to 
find the displacement (domain velocity) field that satisfies the required 
regularity conditions. Thus 

CK]{V} = {f} [8.1] 

where [K] is the reduced stiffness matrix, {V} is the velocity vector of 
the nodes of varying domain, and {f} is the unknown ficticious boundary 
force that produces a perturbation of the boundary. In segmented form, 

Eq. 8.1 becomes 



[ 8 . 2 ] 


where {V^} is the given perturbation of nodes on the boundary, {V^} is the 
node velocity vector in the interior of the domain and {f^} is the 
ficticious boundary force acting on the varying boundary. Equation for 
the unknown interior node velocity vector can be obtained from Eq. 8.2 as 


l K ddi IV - - i K bdi (V 


[8.3] 




If Bezier curves or surfaces [17] are used for boundary representa- 
tion, positions of control points are selected as design parameter , i * 
1,2,. ..,k. To use Eqs. 4.15, 4.22, and 4.26 for sensitivity computation, 
interior node velcity vector {V^} should be expressed in terms of 
variations of design parameter fib^, i = 1,2,. ..,k. To obtain this 
expression, boundary perturbation {V^} should be written in terms of 
variation fib of the design parameter. Once the inverse matrix [K dd ] * is 
obtained, Eq. 8.3 can be used to express {V d J in terms of the variation 
fib of design parameter. However, this requires large computational 
effort. To gain computational efficiency, following method is used: 


first by pertubing a design parameter b^ a unit magnitude, boundary 
perturbation {V^} can be obtained. Then Eq. 8.3 can be solved to 


obtain {V d }. Using {V d } and displacement shape functions, Eqs. 4.15, 
4.22, or 4.26 can be evaluated which gives . This method requires to 
solve Eq. 8.3 k times. However, much as in ti\e adjoint analysis, this is 


an efficient calculation, if Eq. 8.3 has already been solved, requiring 


only evaluation of the solution of the same set of finite element 
equations with different right side for each unit perturbation of b^ , i = 
1,2,.. .,k. Once design change has been determined using iterative design 
process, regridding of interior grid points can be carried out using { V d } 
of Eq. 8.3. 

The automatic regridding method presented here can be used with the 
boundary-layer approach very effectively. That is, for the fixed domain 

V d can be set equal to zero and thus reduce the dimension of [ K dd ] in 
Eq. 8.3. 

To demonstrate feasibility of the method, an engine bearing cap shown 


in Fig. 8.1 is treated. The engine bearing cap is modeled as a three 
dimensional elastic solid. Due to symmetry, only the right half of the 


cap is analyzed. The finite element configuration and loading conditions 
are shown in Fig. 8.1. The material used is steel with Young's modulus 
and Poisson's ratio of E = 1.0 x 10^ psi and v = 0.3, respectively. The 
finite element model shown in Fig. 8.1 contains 82 elements, 768 nodal 


points, and 2111 degrees-of-freedom. For analysis, ANSYS finite element 
STIF95 [26], which is a 20-noded isoparametric element, is used. 

The design variables for this problem are the shape of the varying 
surface r^, distance of clamping bolt center line AB, and distance Cg 
of edge from cap centerline (Fig. 8.2). For surface r^, a Bezier surface 



CLAMPING BOLT FORCE =14,775 lb. 



OIL FILM PRESSURE = 5000 psi 
Figure 8.1 Engine Bearing Cap 



Figure 8.2 Shape Design Parameters of Engine Bearing Cap 

with 4x4 control points is used. For simplicity, only X 2 ~coordi nates of 
four control points Cj thru C 4 are allowed to be varying. That is, 
surface has curvature in the x^-direction only. 

The expression for design sensitivity of von Mises stress averaged 
over individual finite element is obtained from Eq. 4.22. Numerical 


computation of design sensitivity information has been carried out using 
ANSYS finite element code [26] and the computational procedure of Section 
6. For computation of domain velocity vector, Eq. 8.2 is solved using 
ANSYS finite element code. 

Numerical results with a 1 % design change are shown in Table 8.1 for 
randomly selected finite elements. Accuracy of design sensitivity is 
excellent except for elements 5, 22, 36, 56, and 57 where the magnitudes 
of actual change are small. 


Table 8.1. Shape Sensitivity of Engine Bearing Cap 


El. 

No. 

t 1 



<P‘ 

(4>'/Arpxl00)% 

1 

9829.4564 

9727.3229 

-102.1335 

-109.7298 

107.4 

2 

9631.2028 

9565.8641 

-65.3387 

-70.6936 

108.2 

4 

10620.3320 

10648.0140 

27.6820 

25.7748 

93.1 

5 

11444.4800 

11448.0190 

3.5390 

0.4482 

12.7 

7 

13584.5710 

13604.0520 

19.4810 

17.5377 

90.0 

8 

13641.4950 

13674.1020 

32.6070 

31.2587 

95.9 

10 

17933.5910 

17964.5170 

30.9260 

29.8750 

96.6 

11 

18498.6300 

18526.7900 

28.1600 

27.5316 

97.8 

13 

21202.2630 

21222.8600 

20.5970 

20.559 

100.8 

14 

34270.5140 

34294.7650 

24.2510 

23.7614 

98.0 

16 

8367.4820 

8152.6800 

-214.8020 

-230.0584 

107.1 

18 

9686.1116 

9652.6069 

-33.5047 

-38.2462 

114.2 

20 

12670.2480 

12634.3500 

-35.8980 

-38.4216 

107.0 

22 

16248.4050 

16256.2990 

7.8940 

6.3254 

80.1 

24 

30901.3690 

30862.5960 

-38.7730 

-36.4113 

93.9 

26 

7311.4083 

6999.4094 

-311.9989 

-321.7022 

103.1 

27 

6857.7422 

6519.3747 

-338.3675 

-344.1961 

101.7 

29 

7917.4211 

7774.3107 

-143.1104 

-148.4639 

103.7 

30 

7234.2502 

7081.2085 

-153.0417 

-159.7947 

104.4 

32 

9528.8157 

9395.4689 

-133.3468 

-137.7809 

103.3 

33 

8521.5076 

8387.1466 

-134.3610 

-141.4600 

105.3 

35 

13328.4650 

13264.9790 

-63.4860 

-59.4243 

93.6 

36 

12091.2310 

12159.1600 

67.9290 

87.2240 

128.4 

38 

24661.4990 

23849.7610 

-811.7380 

-842.1029 

103.7 

39 

44231.0680 

42109.0220 

-2122.0460 

-2222.5504 

104.7 

41 

7349.5330 

6999.9927 

-349.5403 

360.6076 

103.2 

42 

6920.5279 

6531.4016 

-389.1263 

-404.7290 

104.0 

44 

5998.6512 

5844.9335 

-153.7177 

-165.1199 

107.4 

45 

5762.5105 

5630.1473 

-132.3632 

-142.4745 

107.6 

47 

7016.8980 

6905.9260 

-110.9720 

-115.6099 

104.2 

48 

6822.9614 

6736.9477 

-86.0137 

-90.5011 

105.2 

50 

9706.9951 

9449.8267 

-257.1684 

-262.1729 

102.0 

51 

9639.7495 

9411.8005 

-227.9490 

-228.3732 

100.2 

53 

13634.1000 

12964.2560 

-669.8440 

-701.6882 

104.8 




Table 8.1 Continued 


54 

19874.8650 

18390.3600 

-1484.5050 

-1586.2730 

106.9 

56 

6080.3933 

6078.1016 

-2.2917 

-3.6525 

159.4 

57 

6121.4120 

6114.6667 

-6.7453 

-8.1242 

120.4 

59 

5832.8322 

5697.4132 

-135.4190 

-139.1451 

102.8 

60 

6266.5615 

6098.1882 

-168.3733 

-172.9635 

102.7 

62 

7041.7283 

6971.4204 

-70.3079 

-79.6051 

113.2 

63 

8230.6127 

8059.7990 

-170.8137 

-188.3943 

110.3 

65 

4816.1908 

4793.3159 

-22.8749 

-24.1961 

105.8 

66 

4787.5653 

4761.5085 

-26.0568 

-27.6278 

106.0 

68 

3537.2024 

3578.9546 

41.7522 

40.3578 

96.7 

69 

3692.9881 

3725.9600 

32.9719 

31.0383 

94.1 

71 

6541.8233 

6585.9308 

44.1075 

45.1422 

102.4 

72 

6643.3182 

6680.1317 

36.8135 

38.0789 

103.4 

74 

3872.7605 

3898.8240 

26.0635 

25.4267 

97.6 

75 

3820.6962 

3843.9362 

23.2400 

22.5210 

96.9 

77 

3918.9608 

4017.5877 

98.6269 

100.0753 

101.5 

78 

3932.8001 

4024.4881 

91.6880 

92.9458 

101.4 

80 

6240.3854 

6285.3485 

44.9631 

46.3209 

103.0 

81 

6158.4620 

6202.3951 

43.9331 

45.2521 

103.0 

i. 

DESIGN COMPONENT 

METHOD FOR 

BUILT-UP STRUCTURES 



In this section, design sensitivity analysis method for built-up 
structures is presented. Both shape and conventional (sizing) design 
variables for components of built-up structures are considered. For 
conventional design sensitivity analysis, distributed parameter structural 
design sensitivity analysis theory of Ref. 2 is used. 

Consider a built-up structure that is made up of m > 1 structural 
components that are interconnected by kinematic constraints at their 
interfaces. Using the principle of virtual work for built-up structures 
[2], one can obtain the variational formulation of the governing 
equations. 


a u,n (z * z) = f0r 311 z 62 


where 


__ m __ 

a u n (z ’ z) = l a i i( z » z ) 

u, “ i=l uV 


'u,fl 


m 

(Z) = l 

i=l 


l i i^ 

u ,n 


[9.1] 


[9.2] 


[9.3] 


and Z is the space of kinematically admissible displacements [2], which is 

defined as the set of displacement fields that satisfy homogeneous 

boundary conditions and kinematic interface conditions between components. 

In Eqs. 9.2 and 9.3, a . . (z,7) and % . .(7) are energy bilinear and 

u'.sr u\cr 

load linear forms of component i with domain ft 1 . Note from Eqs. 9.2 and 
9.3, the energy bilinear and load linear forms of Eq. 9.1 are simply 
summations of corresponding terms from each component. Thus, as will be 
seen later, the design sensitivity analysis of the built-up structure is a 
simple additive process. 

In this section, design sensitivity information for displacement 
functional is derived for general built-up structures. Once this is done, 
extention to locally averaged stress functional can be carried out easily. 

Define z as the total variation of z, due to both conventional and 
shape design changes [2], 


z ■ -37 z T (x+xV(x), u+t6u) 


T=0 


= ^7 z(x,u+t6u) + ^ z t (x+tV(x),u) 


T=0 


The first variation of Eq. 9.1 is [2] 


[9.4] 


a 6u,fl (z ’ z) + a i,V (z ’ z) + a u,Q (i> ^ = £ 6u,n (z) + *i,V (z) ’ for a11 z e Z 

[9.5] 

where z * [w*, w Z , 6, v]^ for beam/truss component, z = [z*, z Z , z^]^ for, 

1 9 T 

three dimensional elastic solid, and z = [w, v 1 , v z ] for plate/plane 
elastic solid component. The notation of Eq. 9.5 is chosen to clearly 
display which variables are held fixed and which vary. 

Consider a displacement functional that defines the displacement z at 

* p 

nodal point x£ R 

A 

♦ ■ // 6(x-x)z(x) dn r g>6 -j 

r l • j 

n r 

where 6 ( x ) is the Dirac delta measure at the origin. Taking the first 
variation of Eq. 9.6, one obtains [2] 

V - // 6(x-x)z(x) dft 


[9.7] 



Define a variational adjoint equation by replacing z in the term on the 
right of Eq. 9.7 by a virtual displacement "a and equate the result to the 
energy bilinear form evaluated at the adjoint variable X; i.e.. 


a fi (A,A) = // 6(x-x)A(x) dfl, for all xeZ [9.8] 

£i r 

Denote the solution of Eq. 9.8 as X. Since z satisfies kinematic boundary 
and interface conditions [2], Eq. 9.8 can be evaluated at T = z and Eq. 

9.5 can be evaluated at 7 = X, to obtain 


m 


+ ' ■ l [*' i 4(A) - a' . - (z , X) ] 

i = l 6u jJl 1 6U 1 .fi 1 


m 


+ I [i'i i(^) - a* (z,X)] 
i = l u 1 ^ 1 u 1 ,V' 


[9.9] 


where the first term on the right is due to conventional design variation 
and the second term is due to shape design variation. Note that Eq. 9.9 
is valid for general built-up structures that are composed of m > 1 
structural components. For explicit expressions of the second term on the 
right of Eq. 9.9, results of Eqs. 5.1, 5.2, 5.4, and 5.7 - 5.9 can be 
used. For the first term on the right of Eq. 9.9, the interested reader 
is referred to Refs. 2 and 18. 

As seen in this derivation, one can systematically organize design 
sensitivity expressions for built-up structures, using the design 
component method. Moreover, one can develop a modular computer program 
that will carry out numerical integration of terms in Eqs. 5.1, 5.2, 5.4, 
and 5.7 - 5.9, using the same shape functions that are employed by the 
finite element analysis of the original structure. The result will then 
be a general algorithm and numerical method for design sensitivity 
analysis that can be implemented with existing finite element codes as 
shown in Section 6. 

To demonstrate accuracy of the design component method, a truss-beam- 
plate built-up structure is treated in this section. Consider the truss- 
beam-plate built-up structure shown in Fig. 9.1. A distributed vertical 
load f (x) is applied to the plates. The points supported by the trusses 
are at the intersections of two crossing beams nearest to the free edges 
of the structure. No external loads are applied to the truss and beam 
components. The plates and beams are assumed to be welded together. 



Coordinates of intersection points of beams and plates are supposed to be 
in the mid-planes of the plates and neutral axes of the beams. Beam 
components have rectangular cross-sections. 



t 


The design variables for this built-up structure are thickness t'J(x) 

of each plate component, width d 1J (XjJ and height b 1J (x^) of each 

longitudinal beam compoenent, width d 1J (xg) and height b 1 ^ ( x^ ) of each 

transverse beam component, cross-sectional areas h^ (k=l,16) of the four 

4-bar truss components, and positions a. (1=1 ,4) and 0. (j=l,4) of 

* J 

transverse and longitudinal beam components, respectively. The lengths of 
the trusses are fixed, but they may change their ground positions, and the 
outside boundary of the entire structure is fixed; i.e., only the 
locations a. and 0., i,j=l,4, of beams are shape variables. Dimensions of 

* J 

the structure and the numbering and spacing of beams in both directions 
are shown in Fig. 9.1. 

For numerical calculations, conventional and shape design sensitivity 
calculations are carried out seperately. For plate components, 12 degree- 
of-freedom non-conforming rectangular elements [27] are used. For beam 
components, Hermite cubic shape functions are used. The finite element 
model used for design sensitivity analysis is shown in Fig. 9.2. Only one 
quarter of the entire structure is analyzed, due to symmetry. A total of 
484 elements, with 1281 degrees-of-freedom, are used to model the built-up 
structure, including 400 rectangular plate elements, 80 beam elements and 
4 truss elements. 

For numerical data. Young's modulus and Poisson's ratio are 3.0x10^ 

psi and 0.3, respectively. The overall dimensions are x l _2 = 15 in. x 

15 in. At the nominal design, beam components are located 3 in. apart. 

Other dimensions of the built-up structure at the nominal design are; 

uniform thickness t = 0.1 in. for plate components, uniform height h = 0.5 

in. and width d = 0.15 in. for beam components, and length Z = 5.364 in. 

and cross-sectional area h = 0.1 in.^ for truss components. A uniform 

distributed load f = 0.1 lb/in. ^ is applied on the plate components. 

In Table 9.1, design sensitivity accuracy results are given for 

several functionals, with 1% uniform change in all conventional design 

variables except the cross-sectional areas of truss components. Design 

11 22 

sensitivity results for displacements, bending stresses a and a at the 
extreme fiber of longitudinal and transverse beam components, and von 
Mises yield stress 
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g(a) = (a 
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[9.10] 
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Figure 9.2 Finite Element Model of A Truss-Beam-Plate 
Built-up Structure 

at the extreme fiber of plate components are given in Table 9.1. Results 
given in Table 9.1 show good agreement between predictions t' and finite 
differences At except for von Mises yield stresses on plate elements 177, 
358, 380 and 400 which are acceptable but not good. However, note that 
these elements have low von Mises yield stress and At is small, compared 
to others, and may not be accurate. 

For shape design sensitivity calculations, since the built-up 
structure is symmetric with respect to the center C, the locations 
a. and 3,, i,j=l,2, of transverse and longitudinal beams, measured from 

* J 

the center C, are taken as design variables. 




Table 9.1. Conventional Design Sensitivity 

of Truss-Beam-Plate Built-Up Structure 

(a) Displacement 


Node 

No. 

t 1 


Alp 


( ip'/AipxlOO)% 

53 

-2.9755E-04 

-2.8723E-04 

1.0315E-05 

9.5761E-06 

92.8 

95 

-2.5890E-04 

-2.4993E-04 

8.9703E-06 

8.3204E-06 

92.8 

113 

-2.2837E-04 

-2.2048E-04 

7.8934E-06 

7.3284E-06 

92.8 

137 

-4.2898E-04 

-4.1320E-04 

1.5781E-05 

1.61 7 IE— 05 

102.5 

179 

-4.0427E-04 

-3.8949E-04 

1.4788E-05 

1.5042E-05 

101.7 

221 

-3.7551E-04 

-3.6186E-04 

1.3649E-05 

1.3747E-05 

100.7 

268 

-1.2025E-04 

— 1 . 1601E— 04 

4.2443E-06 

4.1932E-06 

98.8 

310 

1.4419E-05 

1.3799E-05 

-6.1962E-07 

-6.1013E-07 

98.5 

335 

1.5392E-04 

1.4824E-04 

-5.6883E-06 

-5.3932E-06 

94.8 

352 

5.5382E-05 

5.3329E-05 

-2.0534E-06 

-1.9995E-06 

97.4 


(b) Bending Stress on Beam Element 

El. 

No. 

* 1 


Aip 


( ip ‘ / A »px 1 00 ) % 

1 

160.415 

155.757 

-4.658 

-5.250 

112.7 

4 

128.072 

124.347 

-3.725 

-4.271 

114.7 

21 

365.447 

355.610 

-9.837 

-9.608 

97.7 

25 

305.855 

297.648 

-8.207 

-8.191 

99.8 

30 

-69.361 

-67.350 

2.011 

2.233 

111.0 

45 

106.547 

103.409 

-3.138 

-3.383 

107.8 

49 

24.570 

23.801 

-0.769 

-0.963 

125.2 

54 

-74.026 

- 72.026 

2.000 

2.141 

107.1 

60 

-4.789 

-4.624 

0.165 

0.162 

98.0 

80 

2.141 

2.063 

-0.079 

-0.080 

102.1 



(c) von Mises 

Stress on Plate Element 


El. 

No. 

'I' 1 

'P 2 

Aip 


U7m*ioo)% 

1 

49.128 

47.808 

-1.320 

-1.560 

118.2 

5 

34.279 

33.323 

-0.956 

-1.059 

110.8 

9 

62.143 

60.507 

-1.636 

-1.695 

103.6 

13 

77.076 

75.067 

-2.009 

-1.982 

98.6 

17 

92.755 

90.404 

-2.351 

-2.157 

91.7 

22 

44.262 

43.040 

-1.222 

-1.459 

119.4 







Table 9.1(c) Continued 


26 

42.663 

41.471 

-1.192 

-1.247 

104.6 

30 

63.522 

61.831 

-1.691 

-1.729 

102.2 

34 

77.930 

75.877 

-2.053 

-1.946 

94.8 

38 

91.868 

89.533 

-2.335 

-2.105 

90.1 

44 

33.332 

32.332 

- 1.000 

-1.189 

118.9 

48 

52.498 

51.074 

-1.424 

-1.477 

103.8 

52 

67.527 

65.707 

-1.820 

-1.761 

96.7 

56 

79.045 

76.974 

-2.071 

-1.895 

91.0 

65 

38.362 

37.273 

-1.090 

-1.150 

105.5 

69 

46.061 

44.733 

-1.328 

-1.297 

97.7 

73 

70.749 

68.874 

-1.875 

-1.784 

95.2 

77 

67.127 

65.280 

-1.847 

-1.609 

87.1 

85 

38.300 

37.201 

-1.099 

-1.104 

100.5 

89 

47.855 

46.506 

-1.348 

-1.288 

95.6 

93 

67.866 

66.083 

-1.783 

-1.759 

98.7 

97 

62.719 

61.070 

-1.649 

-1.602 

97.2 

106 

43.048 

41.825 

-1.223 

-1.191 

97.4 

112 

62.228 

60.640 

-1.588 

-1.477 

93.0 

118 

59.270 

57.835 

-1.435 

-1.385 

96.5 

128 

55.122 

53.691 

-1.431 

-1.362 

95.1 

134 

48.012 

46.800 

-1.211 

-1.119 

92.4 

140 

56.823 

55.539 

-1.284 

-1.214 

94.5 

151 

53.218 

51.887 

-1.331 

-1.183 

88.9 

155 

37.003 

36.148 

-0.855 

-0.735 

85.9 

159 

37.897 

37.097 

-0.800 

-0.691 

86.4 

169 

52.077 

50.777 

-1.300 

-1.163 

89.5 

173 

56.294 

54.890 

-1.404 

-1.292 

92.0 

177 

23.383 

22.944 

-0.439 

-0.276 

62.8 

192 

61.326 

59.720 

-1.605 

-1.568 

97.7 

198 

27.041 

26.479 

-0.562 

-0.451 

80.3 

212 

67.683 

65.843 

-1.840 

-1.854 

100.8 

216 

53.350 

52.077 

-1.273 

-1.264 

99.3 

235 

82.704 

80.664 

-2.041 

-2.127 

104.2 

239 

100.980 

98.872 

-2.107 

-2.104 

99.8 

255 

79.509 

77.476 

-2.033 

-1.968 

96.8 

259 

101.657 

99.539 

-2.117 

-2.040 

96.4 

274 

67.361 

65.559 

-1.802 

-1.765 

98.0 

278 

64.853 

63.399 

-1.454 

-1.403 

96.5 

288 

37.003 

36.148 

-0.855 

-0.735 

85.9 

319 

27.894 

27.305 

-0.590 

-0.586 

99.4 

337 

21.848 

21.459 

-0.389 

-0.419 

107.6 

358 

14.556 

14.427 

-0.129 

-0.190 

146.8 

380 

12.189 

12.077 

-0.113 

-0.161 

142.8 

400 

8.471 

8.381 

-0.090 

-0.119 

132.8 


As mentioned in Section 8, for shape design sensitivity calculations, 
one must define a velocity field that has ^-regularity with its second 
derivative integrable. The beam components are allowed to move in 
transverse directions only. Hence, V 1 is a function of only and V is 
a function of x 2 only. The velocity field in each plate component is 


represented by Hermite cubic functions in each direction. That is, V 1 ( ) 
and V 2 ^) are represented by Hermite cubic functions. To see the 
velocity field representation graphically, consider Fig. 9.3, in which the 
shape functions for V 1 ^) are plotted. In Fig. 9.3, 60 ^ and 6 a 2 denote 
perturbations of locations of transverse beams. From Fig. 9.3, one 
obtains V 1 ^) = <t> 1 (x 1 ) + <|> 2 ( x 1 ) . That is, 

f 2x^ ^ a l 

3 (x ?p-) 60 ^, 0 < x < 

°1 


V 1 (x 1 ) = 


2 (x-a,) Z 3(02-0.) 

] 3 l( x " a l) 1 ] +5^, Oj < X < o 2 


2 L-. 

2 (x-cur r 34-a,) L 

— [(x-a 2 ) 1 ] fidg + a 2 < x < -j 

4 - 

^ [9.11] 


and a similar expression for V 2 (x 2 ). 



Figure 9.3 Shape Functions For The Velocity V*(xj) 

In Table 9.2, design sensitivity accuracy results are given for 
several functionals, with a 0.25% uniform change in shape design 
parameters. Sensitivity results for displacements, bending stress for 
beam components, and von Mises yield stress for plate components, are 
given in Table 9.2. Results given in Table 9.2 show excellent agreement 
between predictions <|>' and the finite differences A«p. 





* 




Table 9.2. Shape Design Sensitivity of 

Truss-Beam-Plate Built-Up Structure 

(a) Displacement 


Node 

No. 

♦ l 

* 2 

Aip 


( ip * / A xpx 1 00 ) % 

53 

-2.9755E-04 

-3.0123E-04 

-3.6843E-06 

-3.6386E-06 

98.8 

95 

-2.5890E-04 

-2.6186E-04 

-2.9640E-06 

-2.9212E-06 

98.6 

113 

-2.2837E-04 

-2.3080E-04 

-2.4214E-06 

-2.3798E-06 

98.3 

137 

-4.2898E-04 

-4.3951E-04 

-1.0535E-05 

-1.0489E-05 

99.6 

179 

-4.0427E-04 

-4.1394E-04 

-9.6675E-06 

-9.6264E-06 

99.6 

221 

-3.7551E-04 

-3.8411E-04 

-8.5986E-06 

-8.5631E-06 

99.6 

268 

-1.2025E-04 

-1.2319E-04 

-2.9391E-06 

-2.9273E-06 

99.6 

310 

1.4419E-05 

1.5210E-05 

7.9096E-07 

7.9419E-07 


335 

1.5392E-04 

1.5978E-04 

5.8589E-06 

5.8686E-06 


352 

5.5382E-05 

5.7951E-05 

2.5690E-06 

2.5786E-06 

100.4 


(b) Bending Stress on Beam Element 

El. 

No. 

t 1 

* 2 

Aip 


(ip'/Aipxl00)% 

1 

160.415 

163.634 

3.219 

3.095 

96.2 

4 

128.072 

131.041 

2.968 

2.871 

96.7 

21 

365.447 

370.311 

4.864 

4.586 

94.3 

25 

305.855 

310.385 

4.530 

4.399 

97.1 

30 

-69.361 

-66.654 

2.708 

2.820 

104.2 

45 

106.547 

109.385 

2.838 

2.794 

98.4 

49 

24.570 

26.757 

2.187 

2.143 

98.0 

54 

-74.026 

-73.243 

0.783 

0.663 

84.7 

60 

-4.789 

-4.853 

-0.064 

-0.065 

101.8 

80 

2.141 

2.113 

-0.028 

-0.026 

91.2 



(c) von Mises 

Stress on Plate Element 


El. 

No. 

* l 




(<|>'/Ai|>xlOO)% 

1 

49.128 

50.043 

0.915 

0.915 

100.1 

5 

31.279 

34.881 

0.602 

0.600 

99.6 

9 

62.143 

63.067 

0.925 

0.924 

99.9 

13 

77.076 

77.949 

0.873 

0.871 

99.8 

17 

92.755 

93.843 

1.088 

1.090 

100.2 

22 

44.262 

45.153 

0.891 

0.891 

100.0 









Table 9.2(c) Continued 


26 

42.663 

43.440 

0.778 

0.776 

99.8 

30 

63.522 

64.437 

0.914 

0.913 

99.8 

34 

77.930 

78.891 

0.961 

0.960 

99.9 

38 

91.868 

93.039 

1.171 

1.172 

100.1 

44 

33.332 

34.143 

0.811 

0.810 

99.9 

48 

52.498 

53.329 

0.831 

0.830 

99.9 

52 

67.527 

68.410 

0.883 

0.881 

99.7 

56 

79.045 

80.149 

1.104 

1.104 

100.0 

65 

38.362 

39.116 

0.754 

0.749 

99.4 

69 

46.061 

46.804 

0.743 

0.742 

99.9 

73 

70.749 

71.663 

0.914 

0.911 

99.7 

77 

67.127 

68.334 

1.207 

1.209 

100.1 

85 

38.300 

39.039 

0.740 

0.738 

99.8 

89 

47.855 

48.520 

0.665 

0.663 

99.7 

93 

67.866 

68.758 

0.892 

0.889 

99.6 

97 

62.719 

64.012 

1.293 

1.294 

100.1 

106 

43.048 

43.773 

0.726 

0.725 

99.9 

112 

62.228 

62.978 

0.749 

0.745 

99.3 

118 

59.270 

60.553 

1.283 

1.284 

100.1 

128 

55.122 

55.881 

0.759 

0.756 

99.7 

134 

48.012 

48.846 

0.834 

0.830 

99.4 

140 

56.823 

58.228 

1.405 

1.408 

100.2 

151 

53.218 

53.793 

0.575 

0.569 

99.0 

155 

37.003 

37.871 

0.868 

0.863 

99.4 

159 

37.897 

39.123 

1.226 

1.226 

100.0 

169 

52.077 

52.764 

0.687 

0.683 

99.3 

173 

56.294 

56.509 

0.215 

0.208 

96.7 

177 

23.383 

24.298 

0.915 

0.907 

99.1 

192 

61.326 

61.258 

-0.068 

-0.075 

109.9 

198 

27.041 

27.324 

0.283 

0.263 

93.0 

212 

67.683 

67.463 

-0.220 

-0.226 

102.7 

216 

53.350 

53.241 

-0.109 

-0.118 

107.7 

235 

82.704 

82.424 

-0.280 

-0.285 

101.6 

239 

100.980 

101.057 

0.077 

0.075 

97.5 

255 

79.509 

79.143 

-0.366 

-0.369 

100.9 

259 

101.657 

101.295 

-0.362 

-0.364 

100.6 

274 

67.361 

66.997 

-0.364 

-0.368 

101.2 

278 

64.853 

64.565 

-0.288 

-0.290 

100.7 

288 

37.003 

37.871 

0.868 

0.863 

99.4 

319 

27.894 

27.717 

-0.177 

-0.182 

102.5 

337 

21.848 

21.623 

-0.225 

-0.230 

102.2 

358 

14.556 

14.247 

-0.309 

-0.313 

101.2 

380 

12.189 

11.934 

-0.255 

-0.255 

100.1 

400 

8.471 

8.315 

-0.156 

-0.154 

99.1 

10. 

AN OPTIMIZATION 

PROBLEM 





To demonstrate application of the design sensitivity analysis method 
for built-up structures in structural design optimization, the truss-beam- 


plate built-up structure of Section 9 is optimized using a sparse matrix 
symbolic factorization technique for iterative structural optimization 
[28] and Pshenichny's linearization method [29]. 

For numerical design sensitivity analysis and optimization, design 
variables are discretized. That is, each plate element has constant 
thickness and each beam element has constant width and height. Since the 
built-up structure is symmetric with respect to the center C, thickness 
t^ , i ■ 1,210, width dj and height b^, i = 1, 40, and the locations c^. , i 
= 1, 2 of transverse and longitudinal beams, measured from the center C, 
are taken as design parameters. Thus, the total number of design 
parameters is 292. 

The optimal design problem of the built-up structure is to minimize 
weight of the structure, subject to the following constraints: 

Displacement at C; ^ = z(C) < 0.105 in. 

Plate element von Mises stress; ^ < 17500 psi , i = 2,211 

Beam element bending stress; -70000 psi < ^ < 70000 psi, i « 212,251 

Plate thickness; 0.05 in. < t- < 0.25 in., i = 1,210 

Beam width; 0.075 in. < d.. < 0.30 in., i = 1, 40 

Beam height; 0.25 in. < b- < 1.00 in., i = 1, 40 

1 L 1 

Beam position; 0 < < — 

Thus, the total number of inequality constraints is 543. 

Same numerical data as in Section 9 is used except weight density is 
0.1 Ib/in.^ and uniformly distributed load f = 17.5 lb/in. ^ is applied to 
the plate components. 

For numerical computation of the shape design sensitivity information, 
derivatives I x , J x , and h x in Eq. 5.1 for beam component and Vt in Eq. 5.8 
for plate component must be computed. Since each plate element has 
constant thickness and each beam element has constant width and height, 
these derivatives are Dirac delta measures and computations of the shape 
design sensitivity information become complicated. 

To avoid this difficulty, the design process is divided into two 
phases. In Phase I, each plate and beam components (not element) have 
constant thickness and constant width and height, respectively. Hence in 
Phase I, the design parameter set includes 6 plate thicknesses, 6 beam 
heights and widths, and 2 beam locations with total of 20 design 
parameters. To assign the same design parameter to elements in a 
component, design variable linking is used in Phase I. Once an optimum 



point is reached in Phase I, the design process is switched to Phase II 
where shape design parameters are fixed and each plate and beam elements 
are allowed to have different design parameters. Hence the total number 
of design parameter is 290 in Phase II. 

For numerical computation, PRIME-750 and Cray-lS computers are used 
for design Phases I and II, respectively. The initial and final designs 
of Phase I is given in Table 10.1. The initial cost is 0.7875 lb and the 
final cost of Phase I is 0.5894 lb. There are 12 design iterations in 

Table 10.1. Initial and Final Designs of Phase I 


Design 


Initial 

Final 

Parameter 

Design 

Design 


t l 


0.0874 

PI ate 

t 2 


0.0500 

Component 

*3 

0.1 

0.0518 

Thickness 

t 4 


0.0501 


H 


0.0815 


t 6 


0.0910 


b l 


0.4234 

Beam 

b 2 


0.6902 

Component 

b 3 

0.5 

0.4990 

Hei ght 

b 4 


0.3974 


b 5 


0.4266 


b 6 


0.4326 


d l 


0.1087 

Beam 

d 2 


0.2214 

Component 

d 3 

0.15 

0.2499 

Width 

d 4 


0.0754 


d 5 


0.0838 


d 6 


0.1286 

Beam 

“l 

1.50 

1.4221 

Position 

a ? 

4.50 

4.5872 



Phase I with average CPU time 15986 seconds per iteration on a PRIME-750 
computer. It is observed in design Phase I that inner beam stiffners move 
inward (c^ = 1.4211 in.) and outer beam stiffners move outward (c^ = 

4.5872 in.). Number of active stress and displacement constraints at the 
final design of Phase I is 31. Thus, it is necessary to calculate 
sensitivity information for 31 constraints out of 251. 

There are 9 design iterations in Phase II with average CPU time 25.84 
seconds per iteration on a Cray-lS computer. Cost function history of 
Phases I and II is shown in Fig. 10.1. The final cost of Phase II is 
0.5388 lb. Number of active stress and displacement constraints at the 
final design of Phase II is 54. A profile of upper half of the final 
design is shown in Fig. 10.2 



Figure 10.1 Cost Function History 




Figure 10.2 A Profile of the Final Design 


V % 
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